The mutual history of Schlegel’s Japanese gecko (Reptilia: Squamata: Gekkonidae) and humans inscribed in genes and ancient literature

Abstract Knowing how the present distribution of organisms was formed is an essential issue in evolutionary ecology. Recently, the distribution of organisms on Earth has been significantly changed by human-mediated dispersal due to globalization. Therefore, significant attention has been paid to such processes. However, although humankind has taken considerable time to achieve modernization, the impact of ancient human activity on ecosystems has not yet been thoroughly studied. We hypothesized that ancient urban development and transitions had a non-negligible effect on species distribution. Inferring the impact of past human activity on ecosystems from ancient literature and verifying that impact by genetic analysis and human history is an effective means of tackling this problem. As geckos, a popular neighbor of human dwellings, are good material for this model, we performed this combination approach using Schlegel’s Japanese gecko, Gekko japonicus. We show that G. japonicus migrated from China to the western Japanese archipelago before Christ. The gecko species dispersed itself from western to eastern the archipelago on a time scale of thousands of years. There are many synchronizations between the dispersal history of G. japonicus and the historical development of human society. It is suggested by such synchronizations that humans have influenced the distribution of G. japonicus many times throughout its dispersal history.


Introduction
How has the current distribution of organisms been affected by human activities? The answer to this question is essential to understanding the origins of biodiversity (1). Sometimes, species have been enabled by human-mediated dispersal to travel much faster and further than their natural dispersal ability (2,3). There are countless organisms whose distributions are artificially expanding worldwide (4,5). However, many of the cases that have received attention are recent events, and there are not many studies from the perspective of evolutionary biology on how human activity has affected organisms over the long history of humans (6)(7)(8).
Schlegel's Japanese gecko, Gekko japonicus (Reptilia: Squamata: Gekkonidae) is widely distributed in eastern China, South Korea, and south of the Tohoku region within the Japanese archipelago. Recently, Japanese populations of this species have been considered to be derived from an old immigrant Chinese population. This hypothesis has been supported by morphological and molecular genetic studies (9,10). It has been suggested by paleospecies distribution modeling that G. japonicus was able to settle in Japan after the middle Holocene when the climate became warmer (10). Due to the habitat preference of G. japonicus for artificial environments, such as buildings and cargo (11), its migration and expansion from China to Japan are predicted to have been affected by ancient human activities.
Here, we conducted a study to clarify the dispersal history of G. japonicus in the Japanese archipelago, thought to have been influenced by human activity over a long-term period. Recently, deciphering historical materials, such as ancient literature, has been emphasized in assessing the effects of humans on past biodiversity and ecosystems (12,13). Gekko japonicus, a popular organism living close to humans, is variously recorded in ancient Japanese literature. Therefore, we first comprehensively deciphered ancient Japanese literature and collected descriptions of reptiles and amphibians (it seems that the two were indistinguishable to ancient peoples). Based on this ancient knowledge, we hypothesized that G. japonicus was introduced to western Japan by at least the 900s AD, and its distribution expanded to eastern Japan relatively recently, after the 1700s AD. Specimens were then collected from various parts of the Japanese archipelago and China, and a population genetic analysis was performed using genotyping by sequencing (GBS). The dispersal scenario estimated from population genetics was compared with the hypothesis based on ancient literature. Finally, the validity of our findings was verified from the viewpoint of human history. Through a combined approach of the humanities and biology, we have clarified the long-term effects of humans on the distribution of G. japonicus.

Scrutiny of ancient literature
Descriptions of lizards, including geckos, have been found in dozens of ancient documents of various ages, and the most significant examples are shown in Fig. 1. As shown in the descriptions in Fig. 1(A), the Japanese name "tokage," which now means "skink," comes from the word for something that is behind a door and may originally have represented a gecko (14,15). The word "tokage" had already appeared in the 10th century (Fig. 1B), suggesting that its existence has been known for a long time (16). We can see from the examples in Fig. 1(C)-(E) that two lizards (skinks and geckos) and one amphibian (newts) were often confused in Japan before 1600 AD, and the notations differed in each document (17)(18)(19)(20)(21). "Geckos" were called "yamori" around the 15th century (22). Since then, the above three have been distinguished, as shown in Fig. 1(F) (23,24). On the other hand, it is suggested from Fig. 1(G) that geckos inhabit areas nearby Kyoto and Kyushu (western Japan), but have not been found in the Kanto region (eastern Japan) as of 1697 (25). This past distribution is shown in Fig. 4(C).

De-novo assembly, phylogenetic analysis, and population structure
The average raw reads were Set. Phylogeny The phylogenetic tree, based on maximum likelihood and the approximate Bayesian method, is shown in Fig. 2(A). Through Model Finder, TVM + F + R10 was selected as the replacement model with the smallest Bayesian information criteria (BIC). In most cases, individuals collected from the same site formed a monophyletic group. Looking at the phylogenetic relationships, the individuals in Nanjing, China, were the most ancestral but did not form an independent monophyletic group. The Kyushu population, excluding Beppu, was the eldest ancestor population of the Japanese archipelago. This clade was regarded as Kyushu 1. Only the Beppu population was phylogenetically far away from Kyushu 1, and that was regarded as Kyushu 2. Still, individuals from Kyushu 1 did also not form a single monophyletic group. Also, while the populations of remote islands, such as Fukue Island and Tsushima, were geographically differentiated from island to island, they were not so differentiated on the mainland of Kyushu. Although we could collect from only one site (Izumo) in the Chugoku region, all the individuals collected from this area formed a monophyletic group. Individuals from Kinki, Kanto, Tokai, and Shikoku formed one large clade. Individuals from the Hokuriku and Tohoku regions formed separate monophyly groups and were geographically differentiated. The support for the nodes in each regional clade, which are color-coded in Fig. 2(A) was Bayes posterior probabilities higher than 0.99 and Ultrafast bootstrap posterior probability (UBPP) (26) higher than 95%.
The results of the ADMIXTURE (27) analysis of the Japanese archipelago are shown in Fig. 2(B). The lowest value of cover error was at K = 5 (0.39297), and the next lowest values were at K = 4, 6, and 7 (0.39369, 0.39637, and 0.39403). Geographically differentiated genetic structures similar to those in the molecular phylogenetic analysis were found. However, at K = 5, the populations of the Shikoku, Kinki, Tokai, and Kanto regions were all grouped into one cluster. This geographical structure is shown in Fig. 2(C). As a result of a detailed ADMIXTURE analysis in the Shikoku, Kinki, Tokai, and Kanto regions, K = 1 had the slightest cover error (0.53364), and the next lowest values, as shown in Fig. 2(D), were at K = 2 (0.53658). The western cluster is the area around the Kofu and Fujinomiya as the eastern end. On the other hand, the western end of the eastern cluster is around Kyoto. In the Tokai region, located between the Kanto and Kyoto regions, it was found that genetic populations were mixed on a gradient.
The observed summary statistics are listed in Table S2, and the posterior distributions for each parameter simulated with ABCtoolbox (28) are in Table S4. A total of two scenarios (M1 and M2) were simulated with ABCtoolbox, and the adopted scenario (M1) is shown in Fig. 3(E). The rejected one (M2) is shown in Figure  S2. The most probable model selected with the multinomial logistic regression (mnlogistic) method was M1 (probability = 0.83). Model M2 (probability = 0.17), in which migration was considered, was rejected. Therefore, migration between regions was considered negligible in subsequent discussions. In model M1, the median age at which each regional group diverged from its ancestral group was as follows: as Nanjing-Fukue Island, 6,506 generations;   Nodes with a closed circle represent Bayes posterior probabilities higher than 0.99 and UBPP higher than 95%. The Fukue Island and Kyushu-T clades are monophyletic groups that were considered separately in the divergence time and demographic history analyses. (B) Genetic structure of the Japanese populations of G. japonicus. Each bar on the horizontal axis represents an individual, and the vertical axis represents the probability that an individual belongs to each cluster (each color). The lowest observed cover error was K = 5. (C) Mapped geographical genetic structure of the Japanese population. (D) Geographical genetic structure of the Shikoku, Kinki, Tokai, and Kanto regions when K = 2. The breakdown of the pie chart represents the probability that an individual belongs to each cluster. The diameter of the pie chart represents the number of individuals at each site.

Genetic structure in the Japanese archipelago
Previous studies based on mitochondrial DNA and microsatellite analyses have suggested that the Japanese populations of G. japonicus have low genetic diversity and poor geographic structure (10). On the other hand, China has more genetic diversity than Japan, even if limited to the east coast. This was explained by the fact that G. japonicus is native to China and exotic to Japan, and gene flow between populations may have been facilitated by the modern development of globalized transportation networks. Since G. japonicus is widely distributed in eastern China, genetic diversity within China may be quite high, but this is not clear at this time. Although it is a matter of concern that only one site in China is used in this study, it is highly likely that the population on the east coast of China, which is geographically closest to the Japanese archipelago, is a lineage close to the ancestral population.
On the other hand, our high-resolution analysis using genomewide SNPs showed geographically differentiated genetic structures of G. japonicus in the Japanese archipelago. In this study, almost all habitats within the Japanese archipelago were covered, and the reliability of the results for the Japanese archipelago is considered to be high. So, the discrepancy with the results of the previous study is simply due to the resolution of the loci analyzed. As distribution tends to be restricted by temperature, some geckos are geographically isolated by barriers, such as mountain ranges, and genetically differentiated even between relatively close land areas (29,30). In the area of the Japanese archipelago considered here, most of the distribution of G. japonicus is in urban areas along the coast. The distribution of the species and its habitat in moun-     tainous regions has not been confirmed (31). Seas and mountains are lined up at the boundaries of each regional group, as shown in Fig. 2(C) and (D). Additionally, as the rejected model M2 of divergence time (Fig. 3F) shows, the degree of gene flow between each was minimal. Therefore, the pattern is similar to the precedents of other gecko species. The divergence of each regional population occurred hundreds to thousands of years ago, suggesting that the influence of gene flow due to the current transportation networks on the expansion of distribution may not be significant. In other words, although genetic exchanges among each regional population have been promoted by modern transport networks, their impact is thought to be minor at the genome-wide level. It is believed that these gene structures were formed by the founder effect, originating from individuals who accidentally human-mediately broke through geographical barriers and were introduced to each region.

Age of migration to the Japanese archipelago and dispersal history
The approximate ages of the bottlenecks of each population were roughly consistent with their divergence times. The demographic history of Aedes albopictus, which is a mosquito that is artificially expanding around the world, is similar to that found in this study. It is known that the individuals of this species experienced a steep bottleneck when they invaded each region (32). Therefore, con-sidering the bottleneck timing as the age of the introduction, the dispersal routes of G. japonicus in the Japanese archipelago are estimated with the divergence scenario ( Fig. 4A and B). It is worth noting that the Tohoku population has experienced two steep bottlenecks in the past. Still, it is thought that the former event reflects the intrusion event into Hokuriku because the individuals from this region are phylogenetically closest to those from the Hokuriku population. Therefore, it is reasonable to think that the actual timing of introduction into Tohoku was the latter. In addition, despite the fact that populations derived from Kyushu in various parts of Japan must have experienced multiple bottlenecks, the dynamics before major bottlenecks were flat in all populations. This can be explained as follows. This estimation was performed based on the genetic information possessed by the current individual. However, much of the genetic information possessed by the ancestral population should have been lost in the most recent major bottleneck. Therefore, it is likely that the dynamics prior to the major bottleneck cannot be plotted in detail. However, the decreasing trend in the starting size of each regional population as one moves to a newer region suggests that the existence of multiple bottlenecks is present.
Based on Japanese human history, we verified whether the human-mediated introduction of G. japonicus in the above scenario was possible. An overview of major Japanese transportation networks from ancient to modern times is shown in Fig. 4(D). The estimated times at which part of the source population in China was introduced to Fukue Island and Kyushu are reasonable. It was suggested that G. japonicus became viable in Kyushu about 8,000 years ago (10). The Goto Islands, including Fukue Island, are a volcanic archipelago, and the last eruption of the volcano on Fukue Island occurred between 2,400 and 2,300 years ago (33). From the results of the demographic history, a pattern was observed in which the population of G. japonicus on Fukue Island immediately returned to its original state after a slight decrease in its effective population size during the approximate time of the eruption. This event, considering this dynamic as a temporary effect of the eruption, is a calibration point that strengthens the validity of the calculation method of the introduction age into each region. Located on the west side of Kyushu, the Goto Islands were gateways to biological exchanges with the continent. It may be that the geckos first invaded the Goto Islands a little earlier than they crossed over to the mainland of Kyushu. This era corresponds to the Jomon period (16,000 to 3,000 years ago) and the Yayoi period (3,000 to 1,700 years ago) in the classification of Japanese human history. Recent archaeological and genetic studies have suggested that rice cultivation had been propagated to Japan from Southeast Asia via China by at least the end of the Jomon period (about 3,000 years ago) (34). It is suggested that there were human exchanges between the Japanese archipelago and the continent during this period, and G. japonicus probably migrated to Kyushu during this time.
The introduction age to the Chugoku region and the Shikoku region corresponds to the middle to late Yayoi period in Japanese human history. Ancient Izumo in the Chugoku region prospered as a major force in Japan at that time, and many archaeological sites have been excavated in the area. Furthermore, ancient Izumo, as the Izumo mythology, appears in the oldest ancient Japanese documents, such as "Kojiki (35)" and "Nihon Shoki (36)." It is said that close exchanges with northern Kyushu through maritime traffic across the Sea of Japan were also related to social development in this region (37)(38)(39). Additionally, maritime traffic across the Seto Inland Sea (the inland sea surrounded by the Kyushu, Shikoku, Chugoku, and Kinki regions) from Kyushu to Kinki has developed since ancient times (40). Gekko japonicus was able to disperse from Kyushu via multiple routes by taking advantage of the ancient development of local communities and cultural exchanges.
This species took a long time to achieve an introduction to the Kinki region, but what kind of insights on that introduction can be obtained from human history? It is said that there was a blank period with almost no records 200 to 400 AD. Still, trade and transportation continued through the Asuka (592 to 710 AD), Nara (710 to 794 AD), and Heian periods (794 to 1185 AD). These were all periods when the capital was located in the Kinki region. It is reasonable to assume, considering the divergence time estimation scenario, that G. japonicus was spread from Shikoku along with the aforementioned maritime traffic in the Seto Inland Sea and the development of human society.
Gekko japonicus seems to have further expanded into the Tokai region after the introduction to the Kinki region had been mostly achieved. A major arterial road called Tokaido was fully developed during this period (41). This was a land route that ran along the coast of the Tokai region and connected the Kinki and Kanto regions. At the end of 1100s AD, the Shogunate (military government) was established in Kamakura, within the Kanto region. This era was the Kamakura period (1185 to 1333 AD). At that time, the Tokaido road functioned as an artery connecting the Shogunate with the capital of the Kinki region. From a topographical point of view, there are few geographical barriers, such as large mountains between the Shikoku, Kinki, and Tokai regions. Consequently, slow natural dispersal was active in addition to humanmediated rapid dispersal. Continuous long-term dispersal events could have created a gradual bottleneck in these regions. By 300 years ago, the geckos had completed their expansion into the Tokai region. Still, individuals of the species seem to have been stopped by the mountains that separate the Tokai and Kanto regions (including Mount Fuji, the highest peak in Japan). It was not until at least 220 years ago, at the end of the Edo period (1603 to 1868 AD), that the gecko was introduced to the Kanto region. This timeline is consistent with the description in "Honcho Shokkan," in which geckos did not inhabit the Kanto region as of 1697 AD. Then, how did they make their way into the Kanto region? Looking at demographic dynamics, the slope of the bottleneck has become steeper since about 150 years ago, during the Meiji period (1868 to 1912 AD). The railway network began to be developed in Japan in the latter half of the 1800s AD. In 1889, the Tokaido Main Line was fully opened, operating a route similar to the Tokaido. Individuals of G. japonicus may have gained the ability to break through the geographical barrier presented by the high mountains by riding this modern pipeline connecting the Tokai and Kanto regions.
The age when individuals of G. japonicus were introduced to the Hokuriku region from the Kinki region by a route different from their spread to the Tokai and Kanto regions was the so-called "Warring States Period" (1467 to 1590 AD) in Japanese history. During that time, wars broke out all over the country. However, this period was also when the movement and distribution of people became active, and the monetary economy developed (42). Possibly, individuals of G. japonicus dispersed in the Hokuriku region with the development of the logistics network between the Middle Ages and the early modern period. Subsequently, these individuals reached Sakata in the Tohoku region. A sea route called "Nishimawari-Kaiun" was established in 1672 AD-merchant ships sailed westward along the coast, stopping by port towns from the Kanto to the Tohoku regions. Sakata, one of the port towns that prospered due to the development of this route (43), functioned as a key route point in the advancement of G. japonicus from the Hokuriku region.

Conclusion
Gekko japonicus migrated from China to Kyushu before prehistoric times, spreading from west to east and south to north in the Japanese archipelago for thousands of years. This becomes clear through an approach in which human history and the natural sciences are integrated. There were many similarities between the spread of G. japonicus and the historical development of Japanese society, suggesting that ancient human-organism interactions are an essential factor in understanding the present distribution of these organisms. The past impact on ecosystems by humans, which has often been overlooked, is more significant than previously imagined. It is hoped that organisms worldwide will be examined from various perspectives to understand the mutual history of humans and nature.

Scrutiny of ancient literature
Using the digital service of the National Diet Library of Japan (https://dl.ndl.go.jp/), in which many originals and manuscripts of ancient Japanese literature are archived, we have scrutinized dictionaries and academic books from various ages, including the most ancient. We comprehensively collected the descriptions of reptiles and amphibians available from these sources.

Sampling
From 2018 to 2021, 182 individuals of G. japonicus were collected from 37 sites in Japan and one in China ( Figure S1). A total of three other species of the genus, Gekko sp. (undescribed species called Nishiyamori in Japanese), G. tawaensis, and G. shibatai, were collected from three other sites in Japan as outgroups. The individual numbers and locations of the samples used in the analysis are summarized in Table S1. To avoid destructive sampling, only autotomized tail tissue was stored in 100% ethanol, and the individuals were released on the spot. It is worth remembering that the tail of the gecko regenerates over time.

Laboratory protocol
Total DNA was extracted from tissue pieces of the gecko tails collected according to the protocol using Nucleospin tissue (TaKaRa, Shiga Pref., Japan). The DNA was further subjected to RNase treatment. Then, double digest restricted-site associated DNA sequencing (ddRAD-seq) was conducted following Peterson et al. (2012) (44). The library was prepared with 40 ng/μL of gDNA for each sample. Samples were digested with EcoRI, MspI, P1, and P2 adapters, and each fragment was ligated. We then pooled them at equimolar concentrations and purified them using the Nucleospin gDNA clean-up kit (TaKaRa). We then selected 300 to 500 base pair (bp) fragments using Pippin Prep (Saga Science, MA, USA). The size-selected DNA fragments were amplified in polymerase chain reactions (PCR) for eight cycles using Phusion PCR reagents (New England Biolabs, MA, USA). Later, we cleaned up the reaction products and removed the PCR primers using the Nucleospin gDNA clean-up kit and Pippin Prep. The tuned ddRAD library was sequenced (150 bp paired-end) using Illumina HiSeqX (Illumina, CA, USA) paired-end sequencing at Macrogen, Japan.

De-novo assembly
Demultiplex was performed using ipyrad 0.9.14 (45) from the sequence, read and divided into individual reads. At this time, a barcode sequence mismatch was not allowed. Then, a section from the 3' end to 5 bp of each lead was trimmed as a quality control if the Qscore was lower than 33. The clustering threshold was set to 85%, the read was clustered for each individual, and a consensus sequence was created. The minimum number of depths allowed was six, and the maximum was 10,000. The ratio of an ambiguous site (N) to a heterozygous base arranged in the consensus sequence was 0.05. The number of raw reads obtained for each individual was unexpectedly large. Therefore, only one of the paired ends (R1) was used as a single end. From the consensus sequence created for each individual, clustering was performed with a threshold of 90% to generate a consensus sequence for the entire population. Finally, a sequence with 180/185 shared loci, a bp number by indel lower than five, and a maximum SNP number per locus lower than 0.2 were extracted from this sequence, and an analysis data set (Set.Phylogeny) was created. In addition, the following datasets were prepared for the other analyses. Set.Structure1 was made from only Japanese individuals, excluding three Chinese and three outgroup individuals with 176/182 shared loci. Set.Structure2 was created from individuals from the Shikoku, Kinki, Tokai, and Kanto regions of Japan, with 73/74 shared loci. Set.Divergence was created by excluding individuals that did not form a monophyletic group in the Kyushu region from Set.Phylogeny with 166/169 shared loci.

Phylogenetic analysis and population structure
A molecular phylogenetic tree based on the maximum likelihood method was created using iqtree 1.6.8 (46) from the loci data of Set.Phylogeny. Model selection was performed based on BIC using Model Finder Plus (47), and the Ultrafast Bootstrap was set to 2,000. We also calculated the posterior probabilities of each node based on the approximate Bayesian method. The genetic structure of the Japanese population was analyzed using ADMIXTURE (27). The data of Set.Structure1 was formatted by PLINK (48) and was used for the analysis. The number of clusters (K) was set to 1 to 30. The data of the Set.Structure2 dataset was used to clarify the more detailed population genetic structure in the Shikoku, Kinki, Tokai, and Kanto regions that formed one cluster in the higherorder structure. Then, the data were analyzed similarly using AD-MIXTURE (the number of K was set to 1 to 15).

Demographic history and divergence time
The demographic history of each regional population was estimated using Stairway plot v2 (49). Stairway plot is a population dynamics estimation software similar to Extended Bayesian Skyline Plot (50). However, it is performed based on the site frequency spectrum (SFS) and is known to provide a more accurate estimation of recent past dynamics than can be made compared to with the latter. First, the observed SFS values for each regional population were projected from the data of Set.Strucuture1 using easySFS (51). Individuals were divided into the following 11 groups based on the phylogenetic and population structure analysis results and geographical divisions (partially different from administrative divisions), such as mountains and the sea. The considered regions were (1)  A preliminary analysis using all individuals in Kyushu as a single population showed almost unchanged population dynamics. Thus, no bottleneck could be detected in Kyushu, even though other regional populations showed results similar to the true analyses shown in Fig. 3(B)-(D). This is the reason why paraphyletic individuals in Kyushu 1 and Kyushu 2 clades with low numbers were excluded, and two monophyletic groups with enough individuals (Kyushu-T and Fukue Island, Fig. 2A) from the Kyushu 1 clade were considered as two regional populations for the true analysis. The number of projections was set to 2N-1. Among the parameters in the input file of Stairway plot, SFS (projected value), nseq (2 N), and nrand were set at unique values for each regional population. Nrand was specified as four ranges (nseq-2)/4, (nseq-2)/2, 3 * (nseq-2)/4, and (nseq-2), as recommended in the manual of Stairway plot software. All other parameters were set to the same values for all regional populations. Following a previous study (29), the mutation rate values were based on those of Cnemaspis geckos (0.025% Myr-1). The exact generation time of G. japonicus is not known, but based on data from a closely related species, G. hokouensis, it is known that G. japonicus individuals take 1 to 2 years to become sexually mature (52). Based on the above, the mutation rate for G. japonicus was converted to 1.5 years per generation (3.75e-8 per site per generation). All other pa-rameters were left as defaults. Analyses were not performed for Kyushu-F, which was removed, or for China, which did not have sufficient populations. Furthermore, divergence time estimation, based on coalescent theory, was performed from the SNP data of the Set.Divergence dataset using ABCtoolbox (28). Fastsimcoal2 (fsc26) (53) was used for the coalescent simulation, and arlsumstat (54) was used to calculate the number of statistical genetic summaries for the simulated sequence. Fst and φ were used as the summary statistics. For the observed data, summary statistics were calculated using arlsumstat (54) from the SNP data of the Set.Divergence dataset. Since the observed data were SNPs, minor allele frequency (MAF) was used as the parameter. As for the other parameters, the number of effective populations and the number of branched generations of each population were used. For each parameter used in the simulation, a noninformation prior distribution was set, and random numbers were generated from the initial distribution for each number of simulations to perform the simulation (Table  S3). In each scenario, 100,000 simulations were performed. Models were selected using the mnlogistic method. The model with the highest probability was adopted. Additionally, the posterior distribution of the parameters was estimated using the rejection method in the adopted scenario. The tolerance rate was set to 0.01 in both cases. The analysis was conducted using the "abc" package (55) in R (56).